##############
## American Identity, Status Threat, and Backlash - Study 1, Public
##############

##############
## Load data
##############
rm(list=ls())
S1_D1 <- read.csv("Study 1 raw.csv", header = TRUE)
head(S1_D1)
nrow(S1_D1)
ncol(S1_D1)


##############
## Exclude participants who failed the attention check
##############
variable.names(S1_D1)
table(S1_D1$Att_Check2)
nrow(S1_D1)

## remove based on main attention check
S1_D2 <- S1_D1[which((S1_D1$Att_Check2) %in% c("meow","Meow","meow meow","meow =^.^=","meow!","meow or purr","Mhew :)","Mew","purr","No noise cat made","meows","meow!","Meow.","no","mew or meow, depending on how silly the cat is.",
                                               "Meow man","pur","MEOW","meow or miaow","meeow","meeowww","purrrrrrr","Meow meow","Meow!","moi","meows","meyow","none","MEOWW","Miya miya","meowwww","MEOW!","meow purrrrrrrrr",
                                               "meow!","mew","mewow","miau"," Meow!","this cat looks irritated, so it would be a growly kind of meow. Give him a treat.",
                                               "PURR", "purrrrr", "Purr","Mrrp?  Mrrrrrow.  mrrrrAAAAAA!","miow","miaow", " Miao or meow","Meow  Prrr Prrr","meow (or purrrrrr...) ","meow (or nya in japanese) ", "meo", "mea", "Meows",
                                               "MEOOOOOW", "Bark", "meow (or nya in japanese)","meow (or purrrrrr...)","Either a purr or a meow", "Meaw", "meow, hiss!","meow, meow, meeeeooww, mew, meow","meow and purr", "Meow Hiss",
                                               "meo","meoww","myaw","meeyow","Miao or meow","Meawww","meowwwwww","Meoww","meow/purr")), ]
nrow(S1_D2)

## check those who failed
S1_D0 <- S1_D1[-which((S1_D1$Att_Check2) %in% c("meow","Meow","meow meow","meow =^.^=","meow!","meow or purr","Mhew :)","Mew","purr","No noise cat made","meows","meow!","Meow.","no","mew or meow, depending on how silly the cat is.",
                                                "Meow man","pur","MEOW","meow or miaow","meeow","meeowww","purrrrrrr","Meow meow","Meow!","moi","meows","meyow","none","MEOWW","Miya miya","meowwww","MEOW!","meow purrrrrrrrr",
                                                "meow!","mew","mewow","miau"," Meow!","this cat looks irritated, so it would be a growly kind of meow. Give him a treat.",
                                                "PURR", "purrrrr", "Purr","Mrrp?  Mrrrrrow.  mrrrrAAAAAA!","miow","miaow", " Miao or meow","Meow  Prrr Prrr","meow (or purrrrrr...) ","meow (or nya in japanese) ", "meo", "mea", "Meows",
                                                "MEOOOOOW", "Bark", "meow (or nya in japanese)","meow (or purrrrrr...)","Either a purr or a meow", "Meaw", "meow, hiss!","meow, meow, meeeeooww, mew, meow","meow and purr", "Meow Hiss",
                                                "meo","meoww","myaw","meeyow","Miao or meow","Meawww","meowwwwww","Meoww","meow/purr")), ]
nrow(S1_D0)
#view(S1_D0$Att_Check2)
library(rstatix)
chisq_test(S1_D0$Att_Check2, S1_D0$Condition) ## inattentiveness does not vary by condition

## remove based on second attention check
table(S1_D2$White_ID_6)
S1_D3 <- subset(S1_D2, S1_D2$White_ID_6==3) 
nrow(S1_D3)
chisq_test(S1_D2$White_ID_6, S1_D2$Condition) ## inattentiveness does not vary by condition

## Exclude all non-whites
table(S1_D3$Ethnicity)
S1_D3 <- subset(S1_D3, S1_D3$Ethnicity==1) 
nrow(S1_D3)

##############
## Demographic summary + balance tests
##############
variable.names(S1_D3)

S1_D3$Condition [S1_D3$Condition == 0] <- "Control"
S1_D3$Condition [S1_D3$Condition == 1] <- "Racial Shift"
S1_D3$Condition [S1_D3$Condition == 2] <- "Common Identity"

is.factor(S1_D3$Condition)

## Demographics
mean(S1_D3$Age, na.rm = TRUE)
sd(S1_D3$Age, na.rm = TRUE)
round(table(S1_D3$Ideology)/nrow(S1_D3)*100,2) ## percentage liberal
round(table(S1_D3$Gender)/nrow(S1_D3)*100,2) ## percentage female

## Balance tests
nrow(S1_D3)
table(S1_D3$Condition)
chisq.test(table(S1_D3$Condition, S1_D3$Age))
chisq.test(table(S1_D3$Condition, S1_D3$Gender))
chisq.test(table(S1_D3$Condition, S1_D3$Education))
chisq.test(table(S1_D3$Condition, S1_D3$Ideology))


##############
## Recode reverse-coded items and rename variables
##############
## recode the second and fourth status threat item
library(car)
S1_D3$Status_Threat_2 = car::recode(S1_D3$Status_Threat_2, "1=7; 2=6; 3=5; 4=4; 5=3; 6=2; 7=1; else = NA")
S1_D3$Status_Threat_4 = car::recode(S1_D3$Status_Threat_4, "1=7; 2=6; 3=5; 4=4; 5=3; 6=2; 7=1; else = NA")

table(S1_D3$American_ID_1)
S1_D3$American_ID_1 = car::recode(S1_D3$American_ID_1, "32=1; 33=2; 34=3; 35=4;36=5; 37=6; 38=7; else = NA")
S1_D3$American_ID_2 = car::recode(S1_D3$American_ID_2, "32=1; 33=2; 34=3; 35=4;36=5; 37=6; 38=7; else = NA")
S1_D3$American_ID_3 = car::recode(S1_D3$American_ID_3, "32=1; 33=2; 34=3; 35=4;36=5; 37=6; 38=7; else = NA")
S1_D3$American_ID_4 = car::recode(S1_D3$American_ID_4, "32=1; 33=2; 34=3; 35=4;36=5; 37=6; 38=7; else = NA")

##############
## Factorial structure outgroup thermometers
##############
S1_D3_thermometers <- S1_D3[,c(22:26)]
head(S1_D3_thermometers)

## Visualize with scree plot
library(ggplot2)
library(psych)
fafitfree <- fa(S1_D3_thermometers,nfactors = 1, rotate = "promax", fm = "ml")
n_factors <- length(fafitfree$e.values)
scree     <- data.frame(
  Factor_n =  as.factor(1:n_factors), 
  Eigenvalue = fafitfree$e.values)
ggplot(scree, aes(x = Factor_n, y = Eigenvalue, group = 1)) + 
  geom_point() + geom_line() +
  xlab("Number of factors") +
  ylab("Initial eigenvalue") +
  labs( title = "Scree Plot", 
        subtitle = "(Based on the unreduced correlation matrix)")


##############
## Reliabilities
##############
colnames(S1_D3)
which( colnames(S1_D3)=="Status_Threat_5" )

psych::alpha(select(S1_D3, 12:16)) ## Status threat
psych::alpha(select(S1_D3, 37:40)) ## American ID
psych::alpha(select(S1_D3, 41:45)) ## White ID
psych::alpha(select(S1_D3, 22:26)) ## Outgroup thermometers


##############
## Mean Scales
##############
S1_D3$mean_Status_Threat <- (S1_D3$Status_Threat_1 + S1_D3$Status_Threat_2 + S1_D3$Status_Threat_3 + S1_D3$Status_Threat_4 + S1_D3$Status_Threat_5)/5
S1_D3$mean_American_Identity <- (S1_D3$American_ID_1 + S1_D3$American_ID_2 + S1_D3$American_ID_3 + S1_D3$American_ID_4)/4
S1_D3$mean_White_Identity <- (S1_D3$White_ID_1 + S1_D3$White_ID_2 + S1_D3$White_ID_3 + S1_D3$White_ID_4 + S1_D3$White_ID_5)/5
S1_D3$mean_Outgroup_Warmth <- (S1_D3$Thermometers_6 + S1_D3$Thermometers_7 + S1_D3$Thermometers_8 + S1_D3$Thermometers_9 + S1_D3$Thermometers_10)/5


##############
## Rescaling
##############
S1_D3$mean_Status_Threat_scaled <- (S1_D3$mean_Status_Threat - min(S1_D3$mean_Status_Threat, na.rm = TRUE)) / (max(S1_D3$mean_Status_Threat, na.rm = TRUE) - min(S1_D3$mean_Status_Threat, na.rm = TRUE))
S1_D3$mean_American_Identity_scaled <- (S1_D3$mean_American_Identity - min(S1_D3$mean_American_Identity, na.rm = TRUE)) / (max(S1_D3$mean_American_Identity, na.rm = TRUE) - min(S1_D3$mean_American_Identity, na.rm = TRUE))
S1_D3$mean_White_Identity_scaled <- (S1_D3$mean_White_Identity - min(S1_D3$mean_White_Identity, na.rm = TRUE)) / (max(S1_D3$mean_White_Identity, na.rm = TRUE) - min(S1_D3$mean_White_Identity, na.rm = TRUE))
S1_D3$mean_Outgroup_Warmth_scaled <- (S1_D3$mean_Outgroup_Warmth - min(S1_D3$mean_Outgroup_Warmth, na.rm = TRUE)) / (max(S1_D3$mean_Outgroup_Warmth, na.rm = TRUE) - min(S1_D3$mean_Outgroup_Warmth, na.rm = TRUE))
S1_D3$IOS_scaled <- (S1_D3$IOS - min(S1_D3$IOS, na.rm = TRUE)) / (max(S1_D3$IOS, na.rm = TRUE) - min(S1_D3$IOS, na.rm = TRUE))
S1_D3$Dems_scaled <- (S1_D3$Thermometers_3 - min(S1_D3$Thermometers_3, na.rm = TRUE)) / (max(S1_D3$Thermometers_3, na.rm = TRUE) - min(S1_D3$Thermometers_3, na.rm = TRUE))
S1_D3$Reps_scaled <- (S1_D3$Thermometers_4 - min(S1_D3$Thermometers_4, na.rm = TRUE)) / (max(S1_D3$Thermometers_4, na.rm = TRUE) - min(S1_D3$Thermometers_4, na.rm = TRUE))
S1_D3$Affirmative_action_scaled <- (S1_D3$Affirmative_action - min(S1_D3$Affirmative_action, na.rm = TRUE)) / (max(S1_D3$Affirmative_action, na.rm = TRUE) - min(S1_D3$Affirmative_action, na.rm = TRUE))

S1_D3$Blacks_scaled <- (S1_D3$Thermometers_6 - min(S1_D3$Thermometers_6, na.rm = TRUE)) / (max(S1_D3$Thermometers_6, na.rm = TRUE) - min(S1_D3$Thermometers_6, na.rm = TRUE))
S1_D3$Asians_scaled <- (S1_D3$Thermometers_7 - min(S1_D3$Thermometers_7, na.rm = TRUE)) / (max(S1_D3$Thermometers_7, na.rm = TRUE) - min(S1_D3$Thermometers_7, na.rm = TRUE))
S1_D3$Hispanics_scaled <- (S1_D3$Thermometers_8 - min(S1_D3$Thermometers_8, na.rm = TRUE)) / (max(S1_D3$Thermometers_8, na.rm = TRUE) - min(S1_D3$Thermometers_8, na.rm = TRUE))
S1_D3$Illegal_scaled <- (S1_D3$Thermometers_9 - min(S1_D3$Thermometers_9, na.rm = TRUE)) / (max(S1_D3$Thermometers_9, na.rm = TRUE) - min(S1_D3$Thermometers_9, na.rm = TRUE))
S1_D3$Legal_scaled <- (S1_D3$Thermometers_10 - min(S1_D3$Thermometers_10, na.rm = TRUE)) / (max(S1_D3$Thermometers_10, na.rm = TRUE) - min(S1_D3$Thermometers_10, na.rm = TRUE))

S1_D3$Ideology_scaled <- (S1_D3$Ideology - min(S1_D3$Ideology, na.rm = TRUE)) / (max(S1_D3$Ideology, na.rm = TRUE) - min(S1_D3$Ideology, na.rm = TRUE))

##############
## Descriptives
##############
library(dplyr)
table(S1_D3$Condition)
S1_D3$Condition <- factor(S1_D3$Condition, levels = c("Control", "Racial Shift", "Common Identity"))
S1_D3$Condition_num <- as.numeric(as.factor(S1_D3$Condition))
table(S1_D3$Condition_num)

S1_t1 <- t(round(S1_D3 %>%
                group_by(Condition_num) %>%
                summarise_at(vars(c(mean_Status_Threat_scaled, mean_American_Identity_scaled, mean_White_Identity_scaled,
                                    mean_Outgroup_Warmth_scaled, IOS_scaled,
                                    Dems_scaled,Reps_scaled,Affirmative_action_scaled,
                                    Blacks_scaled,Asians_scaled,Hispanics_scaled,Illegal_scaled,Legal_scaled)), list(mean = mean, sd = sd), na.rm = TRUE),2))
write.table(S1_t1, file = "S1_t1.txt", sep = ",", quote = FALSE, row.names = F)




##############
## Simple mean differences (descriptive, not used for hypothesis testing)
##############
## Status threat
S1_aov1 <- aov(mean_Status_Threat_scaled ~ factor(Condition), data = S1_D3)
summary(S1_aov1)
S1_reg1 <- lm(mean_Status_Threat_scaled ~ relevel(factor(Condition), ref = 3), data = S1_D3)
summary(S1_reg1)

library(effectsize)
cohens_f(S1_reg1, partial = FALSE)

## American identity
S1_aov2 <- aov(mean_American_Identity_scaled ~ factor(Condition), data = S1_D3)
summary(S1_aov2)
S1_reg2 <- lm(mean_American_Identity_scaled ~ relevel(factor(Condition), ref = 3), data = S1_D3)
summary(S1_reg2)
cohens_f(S1_reg2, partial = FALSE)

## White identity
S1_aov3 <- aov(mean_White_Identity_scaled ~ factor(Condition), data = S1_D3)
summary(S1_aov3)
S1_reg3 <- lm(mean_White_Identity_scaled ~ relevel(factor(Condition), ref = 3), data = S1_D3)
summary(S1_reg3)

## Outgroup warmth
S1_aov4 <- aov(mean_Outgroup_Warmth_scaled ~ factor(Condition), data = S1_D3)
summary(S1_aov4)
S1_reg4 <- lm(mean_Outgroup_Warmth_scaled ~ factor(Condition), data = S1_D3)
summary(S1_reg4)

## Inclusiveness
S1_aov5 <- aov(IOS_scaled ~ factor(Condition), data = S1_D3)
summary(S1_aov5)
S1_reg5 <- lm(IOS_scaled ~ factor(Condition), data = S1_D3)
summary(S1_reg5)

## Democrat warmth
S1_aov6 <- aov(Dems_scaled ~ factor(Condition), data = S1_D3)
summary(S1_aov6)
S1_reg6 <- lm(Dems_scaled ~ factor(Condition), data = S1_D3)
summary(S1_reg6)

## Republican warmth
S1_aov7 <- aov(Reps_scaled ~ factor(Condition), data = S1_D3)
summary(S1_aov7)
S1_reg7 <- lm(Reps_scaled ~ factor(Condition), data = S1_D3)
summary(S1_reg7)

## Affirmative action
S1_aov8 <- aov(Affirmative_action_scaled ~ factor(Condition), data = S1_D3)
summary(S1_aov8)
S1_reg8 <- lm(Affirmative_action_scaled ~ factor(Condition), data = S1_D3)
summary(S1_reg8)

## Separate outgroups
S1_aov9a <- aov(Blacks_scaled ~ factor(Condition), data = S1_D3)
S1_aov9b <- aov(Asians_scaled ~ factor(Condition), data = S1_D3)
S1_aov9c <- aov(Hispanics_scaled ~ factor(Condition), data = S1_D3)
S1_aov9d <- aov(Illegal_scaled ~ factor(Condition), data = S1_D3)
S1_aov9e <- aov(Legal_scaled ~ factor(Condition), data = S1_D3)
summary(S1_aov9a)
summary(S1_aov9b)
summary(S1_aov9c)
summary(S1_aov9d)
summary(S1_aov9e)

##############
## Moderating ideology
##############
S1_aov1_ideo <- aov(mean_Status_Threat_scaled ~ factor(Condition)*Ideology_scaled, data = S1_D3)
summary(S1_aov1_ideo)
S1_reg1_ideo <- lm(mean_Status_Threat_scaled ~ factor(Condition)*Ideology_scaled, data = S1_D3)
summary(S1_reg1_ideo)

S1_aov2_ideo <- aov(mean_American_Identity_scaled ~ factor(Condition)*Ideology_scaled, data = S1_D3)
summary(S1_aov2_ideo)
S1_reg2_ideo <- lm(mean_American_Identity_scaled ~ factor(Condition)*Ideology_scaled, data = S1_D3)
summary(S1_reg2_ideo)

S1_aov3_ideo <- aov(mean_Status_Threat_scaled ~ mean_American_Identity_scaled*Ideology_scaled, data = S1_D3)
summary(S1_aov3_ideo)
S1_reg3_ideo <- lm(mean_Status_Threat_scaled ~ mean_American_Identity_scaled*Ideology_scaled, data = S1_D3)
summary(S1_reg3_ideo)

S1_reg4_ideo <- lm(mean_White_Identity_scaled ~ factor(Condition)*Ideology_scaled, data = S1_D3)
summary(S1_reg4_ideo)

##############
## Mediation through status threat
##############
library(lavaan)
## Dummy-code conditions
S1_D3$Condition_Control <- ifelse(S1_D3$Condition == "Control", 1, 0)
S1_D3$Condition_Shift <- ifelse (S1_D3$Condition =="Racial Shift", 1, 0)
S1_D3$Condition_CIIM <- ifelse(S1_D3$Condition == "Common Identity", 1, 0)

## Comparison shift and CIIM vs. control
## insert DV as needed
S1_D3_simple_mediation <- '
## Direct effect condition on outcomes (path c) 
mean_Outgroup_Warmth_scaled ~ c1*Condition_Shift + c2*Condition_CIIM

## Indirect effect 1: condition on reported American identity (path a)
mean_Status_Threat_scaled ~ a1*Condition_Shift + a2*Condition_CIIM

## Indirect effect 2: reported American identity on status threat + controls (path b)
mean_Outgroup_Warmth_scaled ~ b1*mean_Status_Threat_scaled

## direct effect
direct_shift := c1
direct_CIIM := c2

## indirect effect (a*b): Sobel test (Delta method)
indirect_shift := a1*b1
indirect_CIIM := a2*b1

## total effect
total_shift := c1 + (a1*b1)
total_CIIM := c2 + (a2*b1)
'

## Comparison shift and CIIM vs. control
S1_D3_simple_mediation_Outgroup_fit <- sem(S1_D3_simple_mediation, se = "bootstrap", bootstrap = 5000, data = S1_D3)
S1_D3_simple_mediation_Inclusiveness_fit <- sem(S1_D3_simple_mediation, se = "bootstrap", bootstrap = 5000, data = S1_D3)
S1_D3_simple_mediation_Democrats_fit <- sem(S1_D3_simple_mediation, se = "bootstrap", bootstrap = 5000, data = S1_D3)
S1_D3_simple_mediation_Republicans_fit <- sem(S1_D3_simple_mediation, se = "bootstrap", bootstrap = 5000, data = S1_D3)
S1_D3_simple_mediation_Affirmative_fit <- sem(S1_D3_simple_mediation, se = "bootstrap", bootstrap = 5000, data = S1_D3)
summary(S1_D3_simple_mediation_Outgroup_fit, fit.measures = TRUE, rsquare = TRUE)
summary(S1_D3_simple_mediation_Inclusiveness_fit, fit.measures = TRUE, rsquare = TRUE)
summary(S1_D3_simple_mediation_Democrats_fit, fit.measures = TRUE, rsquare = TRUE)
summary(S1_D3_simple_mediation_Republicans_fit, fit.measures = TRUE, rsquare = TRUE)
summary(S1_D3_simple_mediation_Affirmative_fit, fit.measures = TRUE, rsquare = TRUE)


## Comparison control and CIIM vs. shift
## insert DV as needed
S1_D3_simple_mediation_2 <- '
## Direct effect condition on outcomes (path c) 
Affirmative_action_scaled ~ c1*Condition_Control + c2*Condition_CIIM

## Indirect effect 1: condition on reported American identity (path a)
mean_Status_Threat_scaled ~ a1*Condition_Control + a2*Condition_CIIM

## Indirect effect 2: reported American identity on status threat + controls (path b)
Affirmative_action_scaled   ~ b1*mean_Status_Threat_scaled

## direct effect
direct_control := c1
direct_CIIM := c2

## indirect effect (a*b): Sobel test (Delta method)
indirect_control := a1*b1
indirect_CIIM := a2*b1

## total effect
total_control := c1 + (a1*b1)
total_CIIM := c2 + (a2*b1)
'


## Comparison control and CIIM vs. shift
S1_D3_simple_mediation_Outgroup_fit_2 <- sem(S1_D3_simple_mediation_2, se = "bootstrap", bootstrap = 5000, data = S1_D3)
S1_D3_simple_mediation_Inclusiveness_fit_2 <- sem(S1_D3_simple_mediation_2, se = "bootstrap", bootstrap = 5000, data = S1_D3)
S1_D3_simple_mediation_Democrats_fit_2 <- sem(S1_D3_simple_mediation_2, se = "bootstrap", bootstrap = 5000, data = S1_D3)
S1_D3_simple_mediation_Republicans_fit_2 <- sem(S1_D3_simple_mediation_2, se = "bootstrap", bootstrap = 5000, data = S1_D3)
S1_D3_simple_mediation_Affirmative_fit_2 <- sem(S1_D3_simple_mediation_2, se = "bootstrap", bootstrap = 5000, data = S1_D3)
summary(S1_D3_simple_mediation_Outgroup_fit_2, fit.measures = TRUE, rsquare = TRUE)
summary(S1_D3_simple_mediation_Inclusiveness_fit_2, fit.measures = TRUE, rsquare = TRUE)
summary(S1_D3_simple_mediation_Democrats_fit_2, fit.measures = TRUE, rsquare = TRUE)
summary(S1_D3_simple_mediation_Republicans_fit_2, fit.measures = TRUE, rsquare = TRUE)
summary(S1_D3_simple_mediation_Affirmative_fit_2, fit.measures = TRUE, rsquare = TRUE)


##############
## Mediation through status threat and national identity + controls
##############
library(lavaan)
## Dummy-code conditions
S1_D3$Condition_Control <- ifelse(S1_D3$Condition == "Control", 1, 0)
S1_D3$Condition_Shift <- ifelse (S1_D3$Condition =="Racial Shift", 1, 0)
S1_D3$Condition_CIIM <- ifelse(S1_D3$Condition == "Common Identity", 1, 0)

## Comparison shift and CIIM vs. control
## insert DV as needed
S1_D3_serial_mediation <- '
Affirmative_action_scaled ~ c1*Condition_Shift + c2*Condition_CIIM

## Indirect effect 1: condition on reported American identity and controls (path a)
mean_American_Identity_scaled ~ a1*Condition_Shift + a2*Condition_CIIM
mean_White_Identity_scaled ~ a3*Condition_Shift + a4*Condition_CIIM

## Indirect effect 2: reported American identity on status threat + controls (path b)
mean_Status_Threat_scaled ~ b1*mean_American_Identity_scaled + b2*mean_White_Identity_scaled

## Indirect effect 3: status threat on DV (path d)
Affirmative_action_scaled ~ d1*mean_Status_Threat_scaled

## direct effect
direct_shift := c1
direct_CIIM := c2

## indirect effect (a*b): Sobel test (Delta method)
Serial_shift_AM_ID := a1*(b1)*d1
Serial_CIIM_AM_ID := a2*(b1)*d1
Serial_shift_WH_ID := a3*(b2)*d1
Serial_CIIM_WH_ID := a4*(b2)*d1

## total effect
total_shift := c1 + (a1*(b1+b2)*d1) + (a3*(b1+b2)*d1)
total_CIIM := c2 + (a2*(b1+b2)*d1) + (a4*(b1+b2)*d1)

## covariances
mean_American_Identity_scaled ~~ mean_White_Identity_scaled
'

## Comparison shift and CIIM vs. control
S1_D3_serial_mediation_Outgroup_fit <- sem(S1_D3_serial_mediation, se = "bootstrap", bootstrap = 5000, data = S1_D3)
S1_D3_serial_mediation_Inclusiveness_fit <- sem(S1_D3_serial_mediation, se = "bootstrap", bootstrap = 5000, data = S1_D3)
S1_D3_serial_mediation_Democrats_fit <- sem(S1_D3_serial_mediation, se = "bootstrap", bootstrap = 5000, data = S1_D3)
S1_D3_serial_mediation_Republicans_fit <- sem(S1_D3_serial_mediation, se = "bootstrap", bootstrap = 5000, data = S1_D3)
S1_D3_serial_mediation_Affirmative_fit <- sem(S1_D3_serial_mediation, se = "bootstrap", bootstrap = 5000, data = S1_D3)
summary(S1_D3_serial_mediation_Outgroup_fit, fit.measures = TRUE, rsquare = TRUE)
summary(S1_D3_serial_mediation_Inclusiveness_fit, fit.measures = TRUE, rsquare = TRUE)
summary(S1_D3_serial_mediation_Democrats_fit, fit.measures = TRUE, rsquare = TRUE)
summary(S1_D3_serial_mediation_Republicans_fit, fit.measures = TRUE, rsquare = TRUE)
summary(S1_D3_serial_mediation_Affirmative_fit, fit.measures = TRUE, rsquare = TRUE)

## Comparison control and CIIM vs. shift
## insert DV as needed
S1_D3_serial_mediation_2 <- '
Reps_scaled ~ c1*Condition_Control + c2*Condition_CIIM

## Indirect effect 1: condition on reported American identity and controls (path a)
mean_American_Identity_scaled ~ a1*Condition_Control + a2*Condition_CIIM
mean_White_Identity_scaled ~ a3*Condition_Control + a4*Condition_CIIM

## Indirect effect 2: reported American identity on status threat + controls (path b)
mean_Status_Threat_scaled ~ b1*mean_American_Identity_scaled + b2*mean_White_Identity_scaled

## Indirect effect 3: status threat on DV (path d)
Reps_scaled ~ d1*mean_Status_Threat_scaled

## direct effect
direct_control := c1
direct_CIIM := c2

## indirect effect (a*b): Sobel test (Delta method)
Serial_control_AM_ID := a1*(b1)*d1
Serial_CIIM_AM_ID := a2*(b1)*d1
Serial_control_WH_ID := a3*(b2)*d1
Serial_CIIM_WH_ID := a4*(b2)*d1

## total effect
total_control := c1 + (a1*(b1+b2)*d1) + (a3*(b1+b2)*d1)
total_CIIM := c2 + (a2*(b1+b2)*d1) + (a4*(b1+b2)*d1)

## covariances
mean_American_Identity_scaled ~~ mean_White_Identity_scaled
'

## Comparison control and CIIM vs. shift
S1_D3_serial_mediation_Outgroup_fit_2 <- sem(S1_D3_serial_mediation_2, se = "bootstrap", bootstrap = 5000, data = S1_D3)
S1_D3_serial_mediation_Inclusiveness_fit_2 <- sem(S1_D3_serial_mediation_2, se = "bootstrap", bootstrap = 5000, data = S1_D3)
S1_D3_serial_mediation_Democrats_fit_2 <- sem(S1_D3_serial_mediation_2, se = "bootstrap", bootstrap = 5000, data = S1_D3)
S1_D3_serial_mediation_Republicans_fit_2 <- sem(S1_D3_serial_mediation_2, se = "bootstrap", bootstrap = 5000, data = S1_D3)
S1_D3_serial_mediation_Affirmative_fit_2 <- sem(S1_D3_serial_mediation_2, se = "bootstrap", bootstrap = 5000, data = S1_D3)
summary(S1_D3_serial_mediation_Outgroup_fit_2, fit.measures = TRUE, rsquare = TRUE)
summary(S1_D3_serial_mediation_Inclusiveness_fit_2, fit.measures = TRUE, rsquare = TRUE)
summary(S1_D3_serial_mediation_Democrats_fit_2, fit.measures = TRUE, rsquare = TRUE)
summary(S1_D3_serial_mediation_Republicans_fit_2, fit.measures = TRUE, rsquare = TRUE)
summary(S1_D3_serial_mediation_Affirmative_fit_2, fit.measures = TRUE, rsquare = TRUE)


##############
## Exploratory: Effect on affective polarization
##############
## Compute AP-score
S1_D3$AP_score <- abs((S1_D3$Thermometers_3 - S1_D3$Thermometers_4))
head(S1_D3[1:6,c("Thermometers_3","Thermometers_4","AP_score")])

## Run regressions
S1_aovAP <- aov(AP_score ~ factor(Condition), data = S1_D3)
summary(S1_aovAP)
S1_regAP <- lm(AP_score ~ relevel(factor(Condition), ref = 3), data = S1_D3)
summary(S1_regAP)
cohens_f(S1_regAP, partial = FALSE)


##############
## Inspection text responses to common identity prime
##############
library(wordcloud)
library(RColorBrewer)
library(wordcloud2)
library(tm)

variable.names(S1_D3)
table(S1_D3$Condition)
Open_Text <- S1_D3$Unique
docs_Open_Text <- Corpus(VectorSource(Open_Text)) ## create a corpus

## clean text
docs_Open_Text <- docs_Open_Text %>%
  tm_map(removeNumbers) %>%
  tm_map(removePunctuation) %>%
  tm_map(stripWhitespace)
docs_Open_Text <- tm_map(docs_Open_Text, content_transformer(tolower))
docs_Open_Text <- tm_map(docs_Open_Text, removeWords, stopwords("english"))

## create a document-term matrix
dtm_Open_Text <- TermDocumentMatrix(docs_Open_Text) 
matrix_Open_Text <- as.matrix(dtm_Open_Text) 
words_Open_Text <- sort(rowSums(matrix_Open_Text),decreasing=TRUE) 
words_Open_Text <- data.frame(word = names(words_Open_Text),freq=words_Open_Text)

## generate the word-cloud
set.seed(4321) # for reproducibility 
wd_low <- wordcloud(words = words_Open_Text$word, freq = words_Open_Text$freq, min.freq = 5,
                    max.words=200, random.order=FALSE, rot.per=0.35,
                    colors=brewer.pal(8, "Dark2"))
df_Open_Text <- data.frame(word = names(words_Open_Text),freq=words_Open_Text)



##############
## End of the syntax.






